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We study the effect that uncertainties in the nuclear spin-dependent structure functions have 
in the determination of the dark matter (DM) parameters in a direct detection experiment. We 
show that different nuclear models that describe the spin-dependent structure function of specific 
target nuclei can lead to variations in the reconstructed values of the DM mass and scattering cross- 
section. We propose a parametrization of the spin structure functions that allows us to treat these 
uncertainties as variations of three parameters, with a central value and deviation that depend on 
the specific nucleus. The method is illustrated for germanium and xenon detectors with an exposure 
of 300 kg yr, assuming a hypothetical detection of DM and studying a series of benchmark points 
for the DM properties. We find that the effect of these uncertainties can be similar in amplitude to 
that of astrophysical uncertainties, especially in those cases where the spin-dependent contribution 
to the elastic scattering cross-section is sizable. 



I. INTRODUCTION 

Direct searches of dark matter (DM) aim to observe 
this abundant but elusive component of the Universe by 
detecting its recoils off target nuclei of a detector (for 
a recent review, see, e.g. Ref. A large number of 
experiments have been taking data in the last decades 
or are currently under construction with this objective, 
leading to a very exciting present situation. 

In fact, some experiments have claimed potential sig- 
nals that could be compatible with the detection of a 
weakly-interacting massive particle (WIMP). This is the 
case of the DAM A collaboration Q, which observed an 
annual modulation in the recoil rate on a Nal target that 
was later confirmed by the upgraded DAM A/LIBRA de- 
tector Q| . Similarly, the CoGcNT collaboration, with 
a germanium target, reported an irreducible excess in 
their data that could point towards very light WIMPs 
Q and also observed an annual modulation effect 0] 
although the latter is not easy to reconcile with the 
DAM A/LIBRA result. Finally, the CRESST-II exper- 
iment, which uses CaWC>4 as a target, also reported 
an excess Q over the expected background. However, 
these observations are in conflict with the negative re- 
sults obtained in searches by other experimental collabo- 
rations. Experiments such as CDMS-II 0,0], XENON10 
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[a], XENON100 [lflEll, EDELWEISS 0, SIMPLE [1J], 
KIMS [3, and a combination of CDMS and EDEL- 
WEISS data 15[ are in strong tension with the regions 
of the parameter space compatible with WIMP signals 
in DAMA/LIBRA or CoGeNT. Moreover, a reanalysis 
of CDMS data has been performed in order to look for 
annual modulation with negative results [16 1. 

The clastic scattering cross-section of WIMPs off nuclei 
can be separated in two components, spin-independent 
(SI) cTq 7 ' N , and spin-dependent (SD) <Jq D ' N , which orig- 
inate from different terms in the Lagrangian describing 
the interaction of a DM particle with quarks. The SI 
term stems from scalar or vector couplings and its con- 
tribution to the total WIMP-nucleus cross-section scales 
as the nucleon number squared, A 2 , whereas the SD term 
originates from axial- vector couplings and its total con- 
tribution to the cross-section off nuclei is only a function 
of the total nuclear angular momentum and the DM spin. 
Thus, the SI term typically dominates for heavy nuclei. 

Constraints are normally expressed in terms of the SI 
and SD components of the WIMP-nucleon elastic cross- 
section, a SI and a SD , respectively To date, the most 
stringent constraints on a SI are those obtained from 
the XENON100 data [11| that exclude SI cross-sections 
above a SI « 2 x 10 -8 pb for a mass around 50 GeV, 
as well as XENON10 |9| and the low-energy reanalysis 
of CDMS-II fl, which dominate for light WIMPs. Re- 
garding the SD contribution, the leading bounds from 
direct detection experiments have been provided by 
XENON [nl (SD cross-section with neutrons, o- SD ' n ) 
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and COUPP [18| and PICASSO jl9| (SD cross-section 
with protons, a SD ' p ) but indirect detection experiments 



such as SupcrKamiokande |2 ( J and IceCube 21], as well 



as searches for mono-jet [22|, 23 1 and mono-photon plus 



missing energy in Tevatron [24| and LHC [25l427| lead 
to even more compelling constraints on cr SD ' p . Larger 
and more sophisticated direct detection experiments are 
currently under development that will be able to explore 
the DM parameter space with unprecedented sensitivity. 
This is the case, for example, of the SuperCDMS and 
XENON1T collaborations, which aim at the construc- 
tion of 1 Ton scale detectors based on germanium and 
xenon, respectively. 

In the light of this promising experimental situation, 
it seems plausible that the DM can be discovered in the 
near future in direct detection experiments. In such an 
event, the study of the signal rate and spectrum (differ- 
ential rate) can be used to determine some of the DM 
properties, namely its mass, m x , and elastic scattering 
cross-section 
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- 30( . The precision of this reconstruction 
is very sensitive to the characteristics of the detector and 
is affected by uncertainties in the parameters describing 
the DM halo, as well as in the nuclear form factors. As- 
trophysical un certain ties have been widely discussed in 
the literature 
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34j and they are known to introduce 



significant errors in the determination of the mass and 
scattering cross-section of DM. Regarding nuclear uncer- 
tainties, those in the SI form factor have been argued to 
be relatively small 35|. The effect of variations in the 
SD form factors has not been previously addressed and 
constitutes the objective of this work. 

We consider the hypothetical future observation of 
a DM candidate in a direct detection experiment 
and, sampling over the three-dimensional space of 
(m x , a , a ), we investigate how the reconstruction of 
these quantities is affected by nuclear uncertainties in the 
spin-dependent structure function of the target nucleus. 
In order to do so, we propose a description of structure 
functions based on three parameters, which enlarge the 
parameter space sampled, and allow us to incorporate 
uncertainties in a consistent and systematic way. This 
provides a general method, applicable to any detector 
target. We particularize our analysis for the case of a ger- 
manium detector (such as, e.g., SuperCDMS), for which 
we consider the spin-dependent structure functions pro- 
vided by the analysis of various groups 0-111, and for 
xenon detectors (such as, e.g., the future XENON1T), for 
which we use the structure functions derived in Ref. 13911 



and [40|. 

We observe that the effect of nuclear uncertainties in 
SD structure functions can lead to variations in the re- 
constructed DM mass and SD elastic cross-section, the ef- 
fect being more important in those scenarios in which the 
SD term in the WIMP-nucleus cross-section is the main 
contribution to the total detection rate. In such cases 
uncertainties in the spin-dependent structure functions 
are similar in amplitude to those induced by astrophys- 
ical uncertainties in the DM halo parameters, although 
the latter also affect the SI component. 

The paper is organized as follows. In Sccfll] we in- 
troduce the formalism used to compute the recoil event 
rate, emphasizing the role of SD interactions. We con- 
centrate on the case of a germanium detector, introduce 
the models available in the literature that describe the 
spin-dependent structure function and comment on their 
differences. Sec. lIIII describes the generation of the simu- 
lated data for a set of benchmark models, and the imple- 
mentation of the scanning algorithm to probe the phe- 
nomenological parameter space. In Sec. lIVI we show the 
reconstruction of DM parameters for each benchmark 
scenario, using different nuclear models for the SD struc- 
ture function and investigating how this alters the pre- 
dictions for the DM properties. In Sec.fV] we present a 
paramctrization of the SD structure function that allows 
us to systematically account for uncertainties when scan- 
ning over our parameter space, and we apply the method 
to the cases of germanium and xenon detectors. Our 
conclusions are summarised in Scc.lVIl 



II. 



NUCLEAR UNCERTAINTIES IN DIRECT 
DARK MATTER DETECTION 



The differential event rate for the elastic scattering of 
a WIMP with mass m x off a nucleus with mass is 
given by 



dR 

cLEr rriN fn 



Po 



vf(v)^(v,E R )dv, (1) 
dEn 



where po is the local WIMP density and f(v) is the 
WIMP velocity distribution in the detector frame nor- 
malized to unity. The integration over the WIMP veloc- 
ity v is performed from the minimum needed to induce 
a recoil of energy E R , v min = y / m N E R /2fi 2 N , to the es- 
cape velocity, v esc , above which WIMPs arc not bound 
to the Milky Way. The WIMP-nucleus elastic scatter- 
ing cross-section, da/dE R , is expressed as a function of 
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the recoil energy, and /ijv = m^m x j '(mjv + Tn x ) is the 
reduced mass. The total event rate is calculated by in- 
tegrating Eq. ([T]) over all the possible recoil energies in a 
window defined by a threshold energy Et and a maximal 
energy E max , both depending on the experiment 

In general, the WIMP-nucleus cross-section is sepa- 
rated into a SI and a SD contribution, as follows: 



da 



m N 



dE R 2p 



AT' 



SI,N, 
'0 



SD, N t?1 
On r • 



SD 



(Er) 



(2) 



where a. 



SI,N 




and a„ D ' N are the SI and SD WIMP- 



nucleus cross-sections at zero momentum transfer. 
Fsi(Er) and Fsd(E r ) are the SI and SD form factors 
that account for the coherence loss which leads to a sup- 
pression of the event rate for heavy WIMPs or heavy 
nuclei. The differential rate, dR/dE R , depends on the 
recoil energy En through the form factors and the mini- 
mal velocity v min (E R ). 

The total number of recoils, as well as their distribu- 
tion in energy, are affected by uncertainties in the nuclear 
form factors (both SI and SD) and in the parameters de- 
scribing the DM halo (usually referred to as astrophysical 
uncertainties). Determining the impact of these is cru- 
cial to understand the capability of a DM experiment to 
reconstruct the properties of the WIMP. 

The role of astrophysical uncertainties has been widely 
addressed in the literature. They are known to signifi- 
cantly affect the reconstruction of both t he ma ss and 
scattering cross-section of the DM 
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43-461 Since 



the subject of our work is to study the effect of nuclear 
uncertainties from the form factors, we do not include 
astrophysical ones. We therefore consider a fixed model 
for the the DM halo, namely the Standard Halo Model 
with a escape velocity ofjJear. = 544 km s _1 , a central 
velocity Vq 



density po 



: 230 km s" 1 [41 
0.4 GcV cm" 3 



51 1 , and a local dark matter 



A. Uncertainties in the SI form factors 

Regarding SI interactions, the so-called Woods-Saxon 
form factor is the Fourier transform of the nucleon dis- 



1 In order to take into account the energy resolution of the detec- 
tor, the differential rate is convoluted with a Gaussian, whose 
standard deviation is a function of the recoil energy, as done in 
Ref. E3. 



tribution function pa{x), 

Fsi(q) = I e-*« x p A (x)d 3 x, 



(3) 



where q = \/2mNE R is the momentum transfer. The 
Fermi distribution is assumed for the nucleon distribu- 
tion, 

1 



Pa(x) oc 



1 



(4) 



exp[(r - Ra)/o\ ' 

where R A = (1.23A 1 / 3 - 0.6) fm, A is the nucleon num- 
ber and a = 0.5 fm the surface thickness of the nucleus. 
Although other parametrizations can be found in the lit- 
erature, the Wood-Saxon form factor provides a good de- 
scription of the nuclear structure for energies in the range 
between 1-100 keV, typical of WIMP scatterings. It has 
been shown in Ref. 35| that the differences in the SI form 
factors due to small deformations of the nuclei can be 
safely neglected. In fact, we have explicitly checked that 
this is indeed the case when using realistic nuclear den- 
sity profiles obtained from a state-of-the art mean field 
calculation. Thus, throughout this paper we consider the 
form factor in Eq. ([3|) with no associated uncertainty. 

B. Uncertainties in the SD form factors 

On the other hand, the effect of uncertainties in the 
SD form factors has not been addressed in the literature. 
The SD contribution to the WIMP-nucleus differential 
cross-section in Eq. (|2|) can be expanded as a function 
of the WIMP couplings to the matrix elements of the 
axial- vector currents in protons (a p ) and neutrons (&„), 

f da \ 16 G 2 F m N (J + 1) 

\dE R J SD nv 2 J 

(a p (S p )+a n (S n )) 2 F 2 SD {E R ), (5) 



where J is the total spin of the nucleus and (S p ) ((S n )) is 
the proton (neutron) spin averaged over the nucleus. The 
SD form factor Fj D {E R ) = S(E R )/S(0), is commonly 
expressed as a decomposition into isoscalar (an = a p +a n ) 
and isovector (a± = a p — a n ) couplings, 

S(q) = alS o{q) + a Q aiS i(q) + a 1 S'ii(g), (6) 

where q is the momentum transfer. The quantities 
Soo(g), ^11(9) an d Soi(q) are the spin-dependent struc- 
ture functions (SDSFs), and are computed using nu- 
clear physics models, whereas the couplings a p and a n 
(and consequently a,Q and a\) are specific of the parti- 
cle physics model for DM and are computed from the 
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diagrams describing the WIMP-nucleon interaction. In 
order to continue with a model independent approach we 
assume a specific relation between a p and a ra , and con- 
sider the casc^l a p /a n = ±1. Under this assumption, 
Eq. (|5|) reduces to 



da 



dE h 




R/ 3D b4Crpmjv 2 Q („\ . a P _ i 



(7) 



The SDSFs 5*00(9) an d <Sii(<z) can be calculated using 
a shell-model (ShM) description of the atomic nucleus, 
where the nuclear spin properties are obtained by the 
wave functions of a few valence nucleons, those which do 
not cancel out the spin of the nucleus in pairs. In par- 
ticular, Soo(q) and Sn(q) are related to the transverse 
electric and longitudinal projections of the axial current. 
To calculate these quantities in the ShM, the nucleons are 
placed in energy levels according to the exclusion princi- 
ple, assuming a particular interaction between nucleons 
(typically a harmonic oscillator potential) and including 
as many excited states as possible, making this kind of 
calculation very difficult. 

ShM calculations are generally more reliable for heavy 
nuclei than for light ones. The same holds for nuclei 
close to magic numbers, elements featuring closed shells 
being more easily modeled. An example is 19 F, that has 
9 protons and 10 neutrons, thus only one proton above a 
magic number. On the other hand, the nucleus of 73 Ge is 
much more difficult to model since it has 32 protons and 
41 neutrons, the nearest closed shell being the one with 
28 nucleons. In this case, deviations of the real nucleus 
from the ShM should be expected, as well as differences 
in the results when different ShMs are used. In the first 
part of the paper we consider the case of germanium, for 
which the only natural isotope that contributes to the 
SD cross-section is 73 Ge. 

In the case of 73 Ge, various ShM calculations are avail- 
able in the literature. We consider two different, com- 
monly used parametrizations , from Ressel et al 
and Dimitrov et al. 
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and in the choice of the nuclear interaction potential, but 
both reproduce the value of the magnetic momentum of 
73 Ge. The SDSFs in both cases can be expressed as a 
function of the adimensional quantity u, related to the 
momentum transfer as u — (qb) 2 /2, where b is the oscil- 
lator size parameter, b = A 1 / 6 . 

The SDSFs for the R-and D-models are plotted as a 
function of u in Fig. [1] by means of red dashed and solid 
red lines, respectively. The left (right) panel refers to the 
case a p /a n = 1 (a p /a n — —1). The vertical, black dashed 
lines indicate the values of u that correspond to the 
WIMP search window that we use in our analysis, from 
a threshold energy of 10 keV, to an energy of 100 keV 
(as currently done in CDMS-II). The dotted blue lines 
indicate a gaussian approximation (see Eq. ([8} below) . 
Finally, the blue areas represent the regions spanned by 
a family of curves, obtained by a parametrization which 
interpolates between the R- and D-modcls that will be 
introduced in Sec.fVl 

The two SDSFs differ in the zero momentum value 
(the R-modcl being larger for the whole energy range of 
interest for direct detection), and also in the shape at 
large energies. They both start as decreasing power-laws 
at low-energy flattening out as u increases. However the 
transition happens sooner for the R-model (around u = 
0.5) than for the D-model. The slope for the D-modcl is 
also slightly steeper than for the R-model, especially in 
Su{q). As we will see in Sec. [IV] these differences play 
an important role when determining the DM parameters. 

There are finally some nuclei for which ShM compu- 
tations of their form factors are not available. In these 
cases an approximation was introduced in Ref. [5ri| that 
works well in the low momentum transfer regime, but 
fails towards larger values of q, 



S ij (q)=S(0)e- 



(8) 



where R, is an effective radius, measured in fm, which 
can be written as, 



[3J, to which we refer as R- and R = 92 A i/3 + 2 68 _ , 78 J( A i/3 _ 3 . 8 )2 + 02 . (9) 



D-models, respectively. They differ in the methodology 



2 This is equivalent to reducing by one the dimensionality of our 
parameter space, assuming a relation between a SD,p and a SD ' n . 
Our analysis can easily be extended to consider the full four- 
dimensional parameter space (m x , a SI , a SD ' p , a SD ' n ), but 
this renders the discussion more cumbersome. Furthermore, par- 
ticle models for DM generally predict |<t SD '"| fs \a SD ' p \. 



III. DETERMINATION OF WIMP 
PROPERTIES 

We consider a set of benchmark scenarios (BM1, BM2 
and BM3) listed in TablcUl that define the phenomeno- 
logical DM parameters (m x , <j SI , a SD ). These bench- 
marks are consistent with possible particle physics mod- 
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u u 

FIG. 1. Spin-dependent structure functions as a function of u, in the case of a v /a„ = 1 (left panel) or a p /a n = — 1 (right 
panel). The solid (dashed) red lines correspond to the D-model (R-model [3]) and the dotted blue line indicates the 
gaussian approximation of Eq. (|8]). The blue region covers the area spanned by the family of curves in Eq. (|14[) . The vertical 
black dashed lines indicate the WIMP search window used in the analysis. 





m x [GeV] a SI [pb] a SD [pb] 


A A s/ X SD 


BM1 
BM2 
BM3 


100 10" 9 10" 5 
50 10" 9 10" 5 
100 10" 9 10" 3 


37.2 36.4 0.8 
42.1 41.2 0.9 
79.6 36.4 43.2 



TABLE I. Phenomenological parameters defining the three 
benchmark models. We include the predicted total number 
of recoil events, A, as well as the number of events (calculated 
using the R-model) A s/ (A SD ) due to SI (SD) interactions, 
for the experimental setup described in the text. 

els for DlvH. We then assume the observation of a DM 
signal in a given direct detection experiment. The differ- 
ential rate is computed for each benchmark point follow- 
ing Eq. (p}, and used to derive the total number of events 
A. 

We first particularize our analysis for the case of a ger- 
manium detector with a total exposure of e — 300 kg yr. 
This could, e.g., correspond to the 1 Ton phase of Super- 
CDMS, operating for a whole year with an efficiency of 
30%. We define the energy window for WIMP searches 
in the range Et = 10 keV and E max = 100 keV, and 
calculate the number of events {A^} in a series of energy 
bins {E l ,E l + AE} with AE = 5 keV. We also include a 
background with a rate of 4 x 10~ 8 days -1 kg -1 keV -1 , 
which is of the order of the background expected for the 



3 In particular, the three benchmarks can be obtained within the 
context of neutralino DM in the general Minimal Supcrsymmct- 
ric Standard Model. 



SuperCDMS experiment in SNOLAB [57| . For the con- 
sidered exposure this means a total of 0.02 background 
events in each of the energy bins considered (i.e., we are 
almost dealing with a background free experiment). We 
assume that this background is flat (energy independent). 

The simulated energy spectra for the three benchmark 
points can be seen in Fig. [2] where the solid red line cor- 
responds to the results when the D-model is used and the 
dashed red line is obtained for the R-model. Practically 
no difference is observed for benchmarks BM1 and BM2. 
This happens because for these two benchmark points 
the main contribution to the total detection rate is due 
to SI interactions. On the contrary, in BM3 the SD con- 
tribution is important and we observe how the predicted 
rate is significantly higher for the R-model than for the 
D-model. This is a consequence of the higher value of Sn 
for the R-model in the whole energy range (see Fig.[T]). 

We treat the quantities {A^} as the experimental in- 
formation from which DM parameters have to be recon- 
structed. Our analysis is based on the Baycs theorem, 
which determines the posterior probability distribution 
(pdf) p(0|D) of a set of parameters (for which a prior 
probability is assumed p(0)) from a set of experimental 
data D, encoded in the likelihood function p(D|0) (or 

£(©)), 



K©|D) 



p(D|0)p(0) 
p(D) 



(10) 



The evidence p(D) in the denominator of Eq. (fTTJ|) is a 
function of only the experimental data. For our purposes 
it works as a normalization factor and can therefore be 
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FIG. 2. Predicted DM spectra for benchmarks BM1, BM2 and BM3 (from left to right) for the experimental setup described 
in the text. The solid and dashed red lines correspond to the predictions using the R- and D-model for the SDSF, respectively. 
The vertical dashed line indicates the energy threshold Et- 



ignored. The pdf in Eq. (fT0|) . in principle, depends on 
the priors p(&) and different choices of priors can affect 
the shape of the final pdf. However, should this hap- 
pen, it would mean that the experimental data are not 
constraining enough and do not dominate the final prob- 
ability distribution. Residual prior dependence can be 
seen, e.g., in Refs. 15814601 . Our scans are performed with 
MultiNest 2.9 6lJ, |62j interfaced with our own code for 
the computation of the number of recoil events and the 
likelihood. Scans are performed with 20000 live points 
and a tolerance of 0.0001. 

In our case, for a given benchmark point the experi- 
mental data is D = ({A;}) and a scan of the parame- 

T SI „SD\ 



ter space = (rn x , a , a ) is performed. The ranges 
considered are m x = 1-10 5 GeV, a SI = 10" 12 -10" 6 pb, 
and a = 10~ 8 — 1 pb. Logarithmic priors are assumed 
for the three variables since the range scanned is quite 
large spanning up to eight orders of magnitude. 

The likelihood £(0) is calculated for each point in the 
scan, computing the number of recoil events Ni in the 
i-th bin, and comparing it with the prediction of the 
benchmark model in the same bin, Ai, assuming that 
each experimental data follows an independent Poisso- 
nian distribution, 



£(©) = n 



N(&) 



A,! 



(11) 



The number of recoil events Ni in the i-th bin are ob- 
tained integrating Eq. ([T]) between Ei and Ei + AE, and 
including a certain number of background events 6; . The 
latter is included as a nuisance parameter in our scans, 
following a Poissonian distribution function with a mean 



value of 0.02. 

The results of our scans are plotted in the next sections 
by means of one- or two-dimensional plots. When the 
probability for a subset of the original is considered, 
one can account for the presence of the hidden parame- 
ters in two different ways: 

• by marginalizing over them, obtaining the pdf for 
the j-th parameter integrating over all the others 



p(9j|D)= f p(®\T>)de 1 ...de j - 1 dG j+1 dQ n ; (12 



• by maximizing over them, obtaining the so-called 
profile likelihood 



6i, 



max 



£(0). 



(13) 



The profile likelihood is usually more sensitive to small 
fine-tuned regions with large likelihood, while the inte- 
gration implemented for the pdf accounts for volume ef- 
fects. Thus, a parameter space characterized by a com- 
plicated likelihood function may result in different pdf 
and profile likelihood distributions for the same param- 
eter. In the following we will present plots for both the 
the pdf and the profile likelihood since they have differ- 
ent statistical meanings and they provide complementary 
information. 



It was recently pointed out in Ref. 34[ that the method 
outlined here is affected by an intrinsic statistical limita- 
tion associated with the fact that only one set of simu- 
lated data is obtained for each benchmark point. In our 
work we do not incorporate this effect, since we want to 
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log(m /GeV) log(rrWGeV) log(m /GeV) 

FIG. 3. Two-dimensional profile likelihood in the plane (m x , <r SI ) for benchmarks BM1, BM2 and BM3, from left to right, 
assuming in the reconstruction that a SD = 0. The inner and outer contours are 68% and 99% confidence level regions, 
respectively. The yellow dot indicates the benchmark values in each case, and the encircled yellow cross marks the positions of 
the best-fit point. 



isolate the variations due to nuclear uncertainties. Fur- 
thermore, our benchmark points are in regions with good 
coverage (see Figs. 2 or 4 in Ref. [3J]) and the number of 
events is relatively large. 

IV. RESULTS 

We begin by considering the three benchmarks BM1, 
BM2 and BM3, with the total rate and energy spectrum 
of nuclear recoils as given in Table[U For concreteness we 
consider explicitly the case a p /a n = — 1, which implies 
looking at only the Sn component of the SDSF, accord- 
ing to Eq.Q. Other choices of a p /a n lead to qualitatively 
similar results (but then a linear combination with the 
other components, Sqo and Soi also appears). We em- 
phasize at this point in that a full analysis can be done 
without fixing this ratio and including a fourth param- 
eter in the scan, however this makes it more difficult to 
isolate the effects of uncertainties in the SDSF. 

As a first exercise in parameter reconstruction, we 
assume (as it is often done) that the SD contribution 
is vanishing and we attempt to reconstruct the WIMP 
mass and SI cross-section from the experimental data (of 
course, this assumption is not made when preparing the 
simulated data from the benchmark points). 

The resulting two-dimensional profile likelihood for 
these quantities are given in Fig. [3] An obvious thing 



to observe is that the reconstruction of these two pa- 
rameters is good for benchmarks BM1 and BM2, since 
in these cases the SI contribution is the dominant term 
in the detection rate. This is obviously not the case in 
BM3, where SD interactions play a more important role. 
For this benchmark point, ignoring the SD contribution 
term leads to an overestimation of the SI independent 
cross-section of approximately a factor two (in order to 
account for the total detected rate). 

Another feature that can be observed, and is consistent 
with the existing literature, is that the goodness of the 
reconstruction is very dependent on the mass of the DM 
(see also Refs. [34|, |4lj, |42(). In partic- 



candidate 28 



ular, we can see how in benchmark BM1 the 99% confi- 
dence level contours arc open for heavy WIMPs, whereas 
this is not the case for BM2 and BM3. In principle, 
increasing the DM mass makes the recoil energy spec- 
trum flatter, as a consequence of the dependence of v m i n 
which enters through the reduced mass. Thus one ex- 
pects to produce a worse fit to the recoil spectrum as we 
scan more massive DM candidates. In benchmarks BM2 
and BM3 this is the reason why heavy masses are dis- 
favoured, however in BM1 the number of events in each 
energy bin is too low to pick up this tendency and very 
massive DM candidates can still produce good fit to the 
data. The presence of a flat background also enhances 
this effect. On the other hand, the contours do not ex- 
tend to low WIMP masses because particles with masses 



4 Note that, since we are scanning over only two parameters, there 



is no need of marginalization or maximization. 
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below m x ^30 — 40 GcV produce a much steeper spec- 
trum. Notice finally that the assumption a SD = leads 
to a lower limit for the SI cross-section that allows us to 
reconstruct the value of a SI up to approximately a factor 
5 (for a fixed value of DM mass). 

An unbiased reconstruction of DM parameters, how- 
ever, has to include the possibility that a SD ^ 0. In 
fact, when we allow for a non-negligible SD contribution 
to the WIMP scattering cross-section, we find that a new 
degeneracy in the parameter space arises: the same de- 
tected rate can be explained by a DM with either pure 
SD or pure SI interactions or, in general, a given com- 
bination of both as we see from Eqs. ([TJ and ©• This 
implies that the closed contours in Fig. [3] can extend to- 
wards arbitrarily small values of <r sl . It is when we ac- 
knowledge this possibility that uncertainties in the SDSF 
play a non-trivial role, as they affect the total rate and 
energy spectrum of WIMP recoils. 

Choosing germanium as a case study, we consider 
the two calculations, R- and D-models, for the SDSF 
of 73 Ge (the isotope that contributes to the SD cross- 
section) that were introduced in Sec. Ill Bl performing the 
Bayesian inference for both SDSFs. Strictly speaking, we 
have to select one SDSF from which the simulated data 
for a given benchmark is generated and, then, one SDSF 
for the computation of the likelihood in Eq. (|TT|) at the 
moment of the scan. The two choices are independent, 
leaving us with four possibilities of combining the two 
SDSFs. In particular, we can either use the same (R- 
or D-) model for the generation of the simulated experi- 
mental data and for the parameter reconstruction or we 
can generate the data with one model and perform the 
reconstruction with the other one. This last possibility 
gives an idea on how sizable the mis-reconstruction of 
DM parameters can be if we use an incorrect model for 
the SDSF, this is, one different to the real one. We per- 
form this computation for the three benchmark points. 

The results for the first benchmark, BM1, are displayed 
in Fig. [4] (profile likelihood) and Fig.[5] (pdf). The four 
rows correspond to the possible combinations of SDSFs 
and we indicate which model is used in the generation 
of data and in the calculation of the likelihood in the 
scan. Comparing the distribution of the profile likeli- 
hood with that of the pdf, we can observe the effect of 
maximization versus marginalization. There are regions 
of the parameter space that are contained in the 99% 
confidence level contour of the profile likelihood which 
are however left out of the credible interval contours of 



the pdf. This happens because the good agreement with 
the data is produced only in a small volume of the three- 
dimensional parameter space and the integration in the 
third dimension decreases the corresponding value for the 
pdf[| 

The plots in the first column of Figs. [4] and [5] repre- 
sent the reconstruction of the WIMP mass and SI cross- 
section, and they can therefore be compared with the 
leftmost plot of Fig. [3] We observe that a SI can now be 
arbitrarily small as long as the a SD is large enough to 
reproduce the observed DM rate and that the assump- 
tion a = leads to over constrained contours. The SD 
form factor (that results for both R-and D-modcls of the 
SDSF) is steeper than the SI form factor. Thus, in prin- 
ciple using the information from the energy spectrum it 
would not be possible to properly fit BM1 data with only 
SD interactions (large values of the WIMP mass provide 
in this case a better fit, since they would correspond to a 
flatter energy spectrum). However, the number of recoils 
events in the high energy bins is too small to be sensitive 
to those differences. 

Similarly, there is no lower bound for the SD cross- 
section. This is evidenced in the second column of both 
Figs. S] and O where a SD is plotted as a function of the 
WIMP mass. As commented above, when trying to fit 
the energy spectrum of a point dominated by SI interac- 
tions (such as BM1) in terms of axial interactions, we get 
better results for large WIMP masses, due to the SDSF 
being steeper. The trend is evident in the pdf plots where 
larger values of masses are associated to a brighter pdf 
than lower masses, but the small number of events pre- 
vents this tendency to have any significant effect on the 
shape of the contours. 

Finally, the degeneracy in the reconstruction of the SI 
and SD contributions for a given set of experimental data 
is clearly evidenced in the third column of Figs.[4]and[5j 
where a SD is plotted versus a and the resulting com- 
patible regions show an "inverse L" pattern. The lower 
left corner of the plot is empty since both the SD and 
SI cross-section are too small to produce the simulated 
number of nuclear recoils, however, as stressed in the 



5 The contours in the two-dimensional plots for the profile like- 
lihood look smoother than for the pdf. This occurs, in particu- 
lar, when integrating over regions with an almost flat likelihood, 
where it is difficult to obtain uniform sampling. The resulting 
pdf can present unphysical structures. The fact that the back- 
ground is included as a free parameter in the scan introduces 
additional fluctuations. 
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FIG. 4. Two-dimensional profile likelihood for the reconstructed parameter space (m x , a , a ) in benchmark model BM1. 
In the first and second rows the model used for the SDSF in the simulated experimental data and in the scan for parameter 
reconstruction is the same (R- model in the first row and D-model in the second). In the third row the D-model is used for the 
simulated data and the R-model for the parameter reconstruction and the reverse is done in the fourth row (see the captions of 
the different panels). The inner and outer contours are 68% and 99% confidence levels, respectively. The yellow dot indicates 
the benchmark value of the parameters, while the yellow encircled cross denotes the position of the best-fit values. 
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FIG. 5. The same as in Fig. [4] but for the marginalized pdf. 
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previous paragraphs, both a s and er can be sizable 
if the WIMP mass is also large. These plots also show 
that for this particular benchmark point SD interactions 
provide the dominant contribution to the WIMP rate for 
a SD > 1(T 3 pb. 

The interpretation of the results for the different rows 
allows us to determine to what extent the uncertainties 
in the SDSF affect the reconstruction of DM parameters. 
We stress again that for BM1 the differential event rate 
is dominated by SI interactions, thus we do not observe 
significant differences when changing the SDSF in the 
computation of the simulated recoils (sec the left panel 
in Fig. [2]). As a consequence, the plots on the first line 
are indistinguishable from those in the third line in Figs.|4] 
and[5l The same happens between the second and fourth 
lines. On the other hand, small differences arise when 
different SDSF are used in the computation of the likeli- 
hood. As already pointed out, the R- and D- models dif- 
fer in the zero-momentum value, as well as in the slope. 
Indeed, we find that when the R-model is used in the 
scan to reconstruct the DM parameters (rows one and 
three) the resulting a SD can be smaller than when the 
D-model is used (rows two and four). This happens be- 
cause the SDSF of the R-model is always larger than in 
the D-model, so the correct number of recoils is repro- 
duced with a slightly smaller a SD . Notice in particular 
how, although the best-fit value for a SD is correctly re- 
constructed when the same SDSF is used for generating 
and reconstructing the points (first two rows), there is a 
mismatch when different models are used. For example, 
if data are generated with the R-model and scanned us- 
ing the D-model (third row), the best-fit value for a SD 
is lower (by about a factor two) than the actual one. Of 
course, the contrary occurs when data are generated with 
the D-model and scanned with the R-model (fourth row) . 
This behaviour can be observed for the three benchmark 
points. Regarding the reconstructed WIMP mass, the 
distribution is similar when either the R- or D-model is 
used, although the latter slightly favours heavier WIMPs 
to compensate for the steeper slope. 

Let us now consider the second benchmark, BM2. We 
proceed as in the previous case and show in Figs.[6]and[7] 
the corresponding reconstruction of the phcnomenologi- 
cal parameters in terms of the profile likelihood and pdf, 
respectively. The difference between profile likelihood 
and pdf (due to the volume effect) is now more strik- 
ing, especially regarding the SD component and WIMP 
mass. As we see in Fig. [6] the regions with a best likeli- 



hood lie around the correct mass but span many orders 
of magnitude in a SD . These regions, however, have a 
small volume and are disfavoured when the pdf is plotted. 
We should emphasize at this point that the information 
from both sources has a different statistical meaning and 
therefore this is no evidence of inconsistency. 

As in the previous scenario, the detection rate in this 
benchmark point is due almost entirely to SI interactions, 
and there are no differences between the simulated data 
with either the R- or D- model for the SDSF (see the 
middle plot in Fig. [2|. However, the number of events is 
now significantly larger and this allows a better determi- 
nation of the slope of the recoil spectrum. This has two 
effects: first, the WIMP mass can be more accurately 
predicted (points with a heavy WIMP being now more 
disfavoured than in the previous example), and second 
leads to larger differences in the reconstruction of a 
when different models for the SDSF are used. Notice, 
for example, how heavy WIMPs are a viable possibility 
only if the contribution from the SD cross-section is suffi- 
ciently large (otherwise the shape of the spectrum is not 
flat enough). In spite of this, the degeneracy between 
a SI and a SD persists. The main difference in the recon- 
struction using the R- or D-model is again the value of 
the lowest <r SD compatible with the data (when a SI is 
negligible), which is smaller for the R- model. Also the 
contours corresponding to the 68% confidence level ex- 
tend towards larger WIMP masses in the case of the re- 
construction using the D-model, in order to compensate 
for its greater steepness. 

Let us finally address the third choice of benchmark 
point, BM3. Contrary to the previous two cases, the 
SD cross-section is a significant contribution to the total 
event rate, as we can see in TableU (although not entirely 
dominant). Therefore we expect that variations in the 
SDSF play a more important role in the reconstruction of 
the DM parameters. The results are displayed in Figs. [5] 
(profile likelihood) and [9] (pdf). 

Ignoring the contribution from the SD term is not a 
good approximation in this scenario. Notice that the re- 
construction of (T in the limit when <r SD = is larger 
(by approximately a factor 2) than the nominal value of 
the BM3 point. This can be appreciated on the lower 
parts of the plots in the third column of both Figs. [5] and 
[9] (and is consistent with the results of Fig. [3]). Moreover, 
the reconstructed value of a also varies, depending on 
whether the scan is performed with the R- or D-model of 
the SDSF, once more due to the different prediction in 




FIG. 6. The same as in Fig.Q]but for benchmark BM2. 




FIG. 7. The same as in Fig.[5]but for the benchmark BM2. 




FIG. 8. The same as in Fig.dbut for benchmark BM3. 
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FIG. 9. The same as in Fig.[5]but for the benchmark BM3. 
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the zero-momentum value. As we already noted in the 
previous benchmark point, an effect in the reconstruction 
of the WIMP mass can also be appreciated between these 
two possibilities. The reconstruction performed with D- 
model favours heavier masses (in fact the 68% confidence 
level contours of the pdf are open for heavy WIMPs) 
than those obtained for R-model, since the D-modcl for 
the SDSF is steeper and this can be compensated with a 
larger value of the WIMP mass, which flattens the spec- 
trum. 

V. PARAMETRIZATION OF UNCERTAINTIES 
IN THE SPIN-DEPENDENT STRUCTURE 
FUNCTIONS 

In the previous section we have shown that the choice 
of model for the SDSF has an important effect in the 
reconstruction of DM parameters. So far our conclusions 
are based on the comparison of the results obtained using 
two different computations for the SDSF of 73 Ge. In 
order to consider these effects in a more systematic way, 
in this section we attempt to include uncertainties in the 
SDSFs as part of the scan. 

To do this, a description of the structure functions has 
to be found in terms of a relatively small number of pa- 
rameters. We propose the use of the following family of 
functions, which reproduces non-trivial features in the 
shape of SDSFs, 

S l3 {u) = N {{1 - ft)e~ au + ft) . (14) 

The parameter N acts as an overall normalization that 
allows us to fit the value at zero-momentum, ft controls 
the height of a possible tail at large momentum and a 
provides the slope of the decreasing part in the the low- 
momentum regim^l 

A. Germanium detectors 

In order to account for uncertainties in the SDSFs we 
have determined the maximum and minimum values of 
the three parameters N , a and ft in Eq. (|14[) which define 



We have explicitly checked that although a five-parameters fit is 
able to reproduce better some features of the SDSF in certain 
nuclei (e.g., 129 Xe and 131 Xe), this has a negligible impact in 
the reconstruction of DM parameters. 



an area that contains the calculations of the R- and D- 
models. The range considered for Sn(q) is the following: 
N = [0.12, 0.21], ft = [0.020, 0.042], and a = [5.0, 6.0]. 
For illustrative purposes we display in Fig.[T] the area 
(in blue) spanned by the family of curves that can be 
obtained by varying the above parameters in the given 
ranges. As we see, the R- and D-models correspond ap- 
proximately to the extremes of the above intervals. 

We repeat the scan for each benchmark extending the 
parameter space including N , a and ft. The number of 
events {A^} of the simulated experimental data are ob- 
tained assuming a SDSF with (N = 0.16, ft = 0.031, a = 
5.5) which is located in the center of the above-mentioned 
ranges. Fig. [10] shows the resulting reconstructed con- 
tours in the profile likelihood of the DM properties in 
the three benchmark models. For comparison, we also 
indicate by means of blue lines the contours of the recon- 
structed DM parameters when nuclear uncertainties are 
not included and where the values of N, a, and ft are 
fixed to their central values. 

We observe that in the case of BM1 the differences 
with respect to the case with no uncertainties are very 
small. One can only observe a slight widening in the de- 
termination of a SD when uncertainties in the SDSF are 
included, but otherwise the reconstructed regions in the 
parameter space show very little differences. This oc- 
curs because in BM1 the DM candidate interacts mainly 
through SI interactions and it is thus fairly independent 
of the details of the SD term. Something similar occurs 
in the case of BM2, although the widening of the re- 
construction of a SD is more evident now. Also the 68% 
confidence level curves corresponding to the WIMP mass 
extend to slightly larger values (notice that the logarith- 
mic scale makes this effect more difficult to observe). Fi- 
nally, it is in benchmark BM3 that the largest effects are 
found, since the SD contribution is larger. Once more, a 
widening in the determination of a SD is observed, which 
is now more evident in the 68% confidence level lines. 
Also the inclusion of uncertainties in the SDSF enlarge 
the contours for large WIMP masses. 

B. Xenon detectors 

The same procedure can be used for xenon detectors. 
Natural xenon contains two isotopes 129 Xe (with a 26.4% 
isotopic abundance) and 131 Xe (21.29%) which are sen- 
sitive to the SD component of the WIMP interaction (in 
particular to the SD cross-section of the WIMP with neu- 
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FIG. 10. Two-dimensional profile likelihood for the reconstructed parameter space (m x , a , a SD ) in benchmark models 
BM1, BM2, and BM3 (from top to bottom), including nuclear uncertainties in the SDSF through the three-parameter model 
introduced in Eq. (|14[) . The inner and outer black contours are 68% and 99% confidence levels, respectively. The solid blue 
line corresponds to the case without uncertainties. The yellow dot indicates the benchmark value of the parameters, while the 
yellow encircled cross the position of the best-fit values. 



trons). As in the case of germanium, we consider various 
parametrizations of the SD form factor for these nuclei 
from Ref. j3i| , in which the nuclear shell model was ap- 
plied to two different potentials describing the nucleon- 



nucleon interaction, the Bonn A 63[ and Nijmegen II 
64J potentials. We also include a recent result from 



Ref. [4(j in which the so called gcn5082 interaction 65 1 
is used. Then we repeat the analysis of the previous 
section modeling the uncertainties in the xenon SDSF 
by means of the parametrizations in Eq. (|14|) , chang- 
ing the values of the (N, a, (3) parameters to define 
the area that contains the above-mentioned models for 
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FIG. 11. The same as in Fig. Q] but for the case of 129 Xe (top row) and 131 Xe (bottom row). The solid (dashed) red lines 
correspond to the ShM calculation using the Bonn A (Nijmegen II) potential [39(]. The solid black line corresponds to the 
determination of Ref. and the dotted black lines are the errors associated to it (the errors for Soo are negligible and are not 
shown) . The dotted blue line indicates the gaussian approximation of Eq. ([HJ . The blue region covers the area spanned by the 
family of curves in Eq. (|14p with the parameters defined in the text. The vertical black dashed lines indicate the WIMP search 
window used in the analysis. 



the SDSFs. In particular, for the S\± component in 
129 Xc we consider N = [0.029, 0.052], a = [4.2, 4.7], 
and (3 = [1.0 x 10~ 3 , 7 x 10" 3 ]. Similarly, in 131 Xe the 
ranges for S n are N = [0.017, 0.027], a = [4.3, 5.0], and 
/3 = [4.2 x 10~ 2 ,6.1 x 10~ 2 ]. The various models for the 
SDSFs are represented in Fig. [11] together with the en- 
velopes for Soo an d Sn in both isotopes. We consider 
the same exposure as in the previous case (e = 300 kg 
yr) but the energy range of the WIMP detection window 
is now taken to be = [8.4, 44.8] kcV, mimicking that 
of the XENON100 experiment. 

Uncertainties in the SDSF for xenon have the same 
qualitative effect as in germanium. Namely, the predic- 
tions for the WIMP mass and the SD component of its 
scattering cross-section are affected. The resulting con- 
tours for the profile likelihood benchmarks BM1, BM2 
and BM3 arc displayed in Figs. [12] We can observe that 



the effect is similar in magnitude to the case of germa- 
nium (despite being a heavier nucleus than germanium, 
the isotopic abundance of the elements sensitive to the 
SD coupling is larger in xenon). Once more, deviations 
are larger for BM2 and BM3 than in BM1. 

The inclusion of uncertainties on SDSF through the 
parametrization in Eq. (fT4|) is a procedure that can be 
applied to other nuclei. In the case of germanium and 
xenon, the existence of different SDSF computations al- 
lowed us to define the ranges in which the three param- 
eters of Eq. (fl4|) are varied. 

C. Comparison with astrophysical uncertainties 

To put our results in context, we need to compare the 
effects of nuclear uncertainties in the SDSF that we just 
discussed with those originating from astrophysical un- 
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FIG. 12. The same as in Fig. [10] but for the case of a xenon detector. 



certainties in the parameters of the DM halo. In order 
to introduce the latter, we have considered a halo model 
motivated by iV-body simulations, which differs from the 
standard halo model in a high- velocity tail H . 6^-6^- 
The distribution function is taken from Ref. |45| and it is 
characterized by the presence of an additional parameter 
k that controls the deviations of F(v) from the standard 



halo model, 
F(v) 



Q(Ve 



■v), 
(15) 

where N k = v$e~ y i fg e dy y 2 { e -^ 2 - y i ) / k - l) fc and 
He = v esc /vQ. In the limit of vanishing k the stan- 
dard halo model is recovered. We then consider that 
the three parameters that define the velocity distribu- 
tion function vary in the ranges v esc G [478, 610] km s _1 , 
V G [170,290] km s _1 , and k e [0.5,3.5], and include 
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them in our scan as nuisance parameters. The local DM 
density is also subject to observational uncertainties. Its 
value can be estimated from a set of experimental con- 
straints that fix the local gravitational potential of the 
Milky Way, with typical values ranging from 0.2 to 0.6 
GeV cm- 3 
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In Figs. [13] and [14] we represent the one-dimensional 
profile likelihood for the DM parameters (m x , a SI , a SD ) 
for benchmarks BM1, BM2 and BM3 in the cases of a 
germanium and xenon detector, respectively. We display 
the reconstruction when no uncertainties are considered 
(blue line), when only nuclear uncertainties in the SDSF 
are included (solid black line) and when only astrophys- 
ical uncertainties are included (dashed black line). As 
noted before, the effect of nuclear uncertainties in the 
SDSF is more evident for BM3 than in BM1 and BM2 
in both germanium and xenon, since in the latter the SD 
component is more important. The prediction for the 
WIMP mass is extended towards larger masses, and as 
we see for BM3 the effect cannot generally be neglected. 
Similarly, the predicted <r SD can vary significantly. In 
BM3 the reconstruction extends towards larger values 
(in BM2 and BM1 the effect is smaller). On the other 
hand, astrophysical uncertainties affect both the recon- 
struction of the three DM parameters, m x , a SI and a 
and are equally relevant, irrespectively of whether the 
main contribution comes from the SD or SI component. 
We can see how nuclear uncertainties generally have a 
smaller effect than astrophysical ones, but they can be 
comparable in some benchmark scenarios, especially re- 
garding the mass reconstruction. This is the case, e.g., 
of BM3 in a xenon detector. 



VI. CONCLUSIONS 

We have studied the effect that uncertainties in the 
nuclear spin-dependent structure functions have in the 
reconstruction of DM properties by means of direct de- 
tection experiments. 

Assuming a hypothetical future observation of DM 
in a direct detection experiment we have systematically 
investigated how well its phenomenological parameters 
(m x , <t si , a SD ) can be determined when uncertainties 
in the SD form factors of the target nuclei are taken into 
account. We focused at first on the case of a germanium 
target and considered two possible models describing the 
SDSF of its isotope 73 Ge, sensitive to SD WIMP cou- 
plings. Using a Bayesian inference algorithm we deter- 



mined for each of these models the pdf and profile likeli- 
hood of the DM parameters in a set of benchmark scenar- 
ios. We observed that if a model is chosen to describe the 
SDSF of a particular nucleus, the reconstruction of the 
DM properties can strongly depend on the choice made 
(see in this sense the comparison between the predictions 
using the R- or D-modcl in Figs.[8]and[9]). In particular, 
differences in the reconstructed values of the WIMP mass 
as well as the SD component of the WIMP-nucleon scat- 
tering cross-section appear. In general these effects are 
more important when the SD contribution to the total 
detection rate is not negligible. 

In the second part of the paper we have proposed a 
description of the SD structure functions in terms of three 
parameters which fit the zero-momentum value and the 
slope of the SDSF, and account for the presence of a high- 
momentum tail. This allows us to include uncertainties 
in the SDSF in the sampling of the parameter space and 
treat them in a consistent and systematic way. Using this 
method we have computed the profile likelihood for the 
DM parameters for the same three benchmark points as 
before, in the case of a germanium-based and a xenon- 
based detector. 

Finally, we have explicitly compared the effect of nu- 
clear uncertainties in the SDSF with those that are as- 
sociated with the parameters of the halo of dark matter. 
We find that uncertainties in the SDSF can even be com- 
parable in magnitude to astrophysical ones when the SD 
contribution to the total detection rate is sizable. 
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FIG. 13. One-dimensional profile likelihood for m x , <r SI , and a SD in BM1, BM2, and BM3 from top to bottom, respectively 
in the case of a germanium detector. The solid blue line corresponds to the case without uncertainties, the black solid line 
represents the results when nuclear uncertainties in the SDSF are included, and the dashed black line denotes the case when 
astrophysical uncertainties are included. The black dot represents the benchmark value of the parameters. 
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FIG. 14. The same as in Fig. 1131 but for the case of a xenon detector. 
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